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ABSTRACT 

The recent implementation of radiative transfer algorithms in numerous hydrodynamics codes has led to 
a dramatic improvement in studies of feedback in various astrophysical environments. However, because of 
methodological limitations and computational expense, the spectra of radiation sources are generally sampled 
at only a few evenly-spaced discrete emission frequencies. Using one-dimensional radiative transfer calcu- 
lations, we investigate the discrepancies in gas properties surrounding model stars and accreting black holes 
that arise solely due to spectral discretization. We find that even in the idealized case of a static and uniform 
density field, commonly used discretization schemes induce errors in the neutral fraction and temperature by 
factors of two to three on average, and by over an order of magnitude in certain column density regimes. The 
consequences are most severe for radiative feedback operating on large scales, dense clumps of gas, and me- 
dia consisting of multiple chemical species. We have developed a method for optimally constructing discrete 
spectra, and show that for two test cases of interest, carefully chosen four-bin spectra can eliminate errors 
associated with frequency resolution to high precision. Applying these findings to a fully three-dimensional 
radiation-hydrodynamic simulation of the early universe, we find that the H II region around a primordial star 
is substantially altered in both size and morphology, corroborating the one-dimensional prediction that discrete 
spectral energy distributions can lead to sizable inaccuracies in the physical properties of a medium, and as a 
result, the subsequent evolution and observable signatures of objects embedded within it. 
Subject headings: dark ages, reionization, first stars — methods: numerical — radiative transfer 



1, INTRODUCTION 

Energy injection by radiative processes fundamentally 
changes the evolution of astrophysical systems, whether it 
be in the context of star formation, galaxy evolution, or 
the growth of super-massive black holes (SMBHs). For in- 
stance, ultraviolet photons from the universe 's first stars (Pop- 
ulation III (PopIII) stars: I Abel et al] 120021) photo-dissociate 
the primary coolant (H2) that first enabled their formation. 
Very recent radiation-hydrodynamic calculations of PopIII 
stars find that PopIII star masses may be limited by proto- 
stellar radiative feedback, perhaps explaining the lack of ev- 
idence fo£_e2COji£_£airinsJability supernovae in the early uni- 
verse (Hosoka wa et al.ll201 ll) . Conventional metal line cool- 
ing driven star formation can be affected by radiative feed- 
back as well. [Krumholz (2006) showed that photo-heating 
around newly formed stars can strongly suppres s fragmen- 
tation in surrounding proto-stellar clouds, while iDale et all 
(2005) see both positive and negative feedback operating 
in radiation-hydrodynamic simulations of star cluster forma- 
tion. Radiative feedback could also be a barri er to efficient 
black hole (BH) growth in the early universe (lAlvarez et al.l 
2009), as X -rays from accreting BHs efficiently photo-heat 
surroundin g gas, leading to smaller Bondi-Hoyle accretion 
rates dBond i & Hovl ilSHl . 

The mere presence of ionizing/dissociating photons ensures 
a change in the chemical and thermal state of a gas, though 
the magnitude of these changes hinges squarely on the num- 
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ber of photons propagating through the gas and their spectral 
energy distribution (SED). Holding the bolometric luminos- 
ity of a radiation source constant, even subtle changes in the 
SED can lead to noticeable differences in the properties of 
the surrounding medium. For example, adjusting the X-ray 
power-law index of a BH accretion spectrum results in ion- 
ization fronts which differ by factors of « 2-3 in radius, and 
temperature p rofiles varying by 10 2 -10 3 K on scales of several 
hundred kpc ( Thom as & Zaroubi 2008). Simply truncating 
the emission of identical X-ray SEDs at harder energies (0.4 
keV rather than 0.2 keV) causes a drastic reduction in heating, 
ionized fraction s, and H2 fractions surro unding 'miniquasars' 
at high redshift dKuhlen & Madaull2005l) . 

Unfortunately, not all radiative transfer algorithms are 
able to represent radiation sources with continuous SEDs, 
or perhaps cannot afford the additional computational ex- 
pense associated with the frequency dependence of the radia- 
tive transfer equation. The natural first step is to represent 
sources as monochromatic emitters, choosing an emission fre- 
quency characteristic of the full SED. Some authors have 
improved upon the monochromatic treatment using 'multi- 
group' methods, which average SED properties and absorp- 
tion cross-sections ov e r one or more frequency bandpasses 
dGnedin& Abelll200U lAubert & Tevssierl 12008b . while oth- 
ers have sampled continuous SEDs at n v frequencies, which 
are generally evenly spaced bins (in linear or log-space) be- 
tween the hydrogen ionization threshold and an upper fre- 
quency cutoff. In either case, there is no clear method of 
deciding how many frequency-averaged bandpasses or dis- 
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crete emission frequencies are required for a given problem, 
and though the standard multi-group treatment is physically 
motivated, it does not guarantee that the photo-ionization and 
photo-heating rates are adequately reproduced as a function 
of column density. 

Frequency resolution has rec ently been s tudied in 
radiation-hydrodynamic se t tings by IWise & A bel (201TJ) and 
IWhalen & Norman] (120081 IWise&Abell (12011b find that for 
the expansion of an H II region around a 10 5 K blackbody 
source in a hydrogen-only medium, the density, temperature, 
velocity, and ionization profiles are well converged for n v > 4. 
Use of a monochromatic spectrum for this problem intro- 
duces significant errors since all photons are absorbed at a 
characteristic column density, whereas multi-frequency treat- 
ments achieve some column density dependent behavior and 
can thus mimic the behavi or of a truly continuous spectrum. 
IWhalen &~N orman (2008) studied the effects of frequency 
resolution in the setting of I-front instabilities, and did not 
achieve convergence until n v > 80 (logarithmically spaced be- 
tween 13.6 and 90 eV). 

The convergence for the test of Wise & Abell (1201 ll) us- 
ing only four frequency bins is reassuring, though the 
prospects for convergence are less clear if one were in- 
terested in the absorption processes of multiple chemical 
species, ionization and heating du e to X-rays and their 
energ e tic secondary photo-electro ns (IShull & van Steenbergl 
119851: iFurlanetto & Stoeverl 120101) . or inhomogeneous me- 
dia. iKramer & Haimanl (I2008L hereafter KH08) briefly com- 
pared monochromatic and continuous treatments of absorbed 
power-law X-ray sources in a study of ionization front thick- 
ness around high-z quasars (the I-front thickness is a po- 
tentially powerful indirect probe of the ionizing spectrum 
of high-z quasars). The hydrogen and helium I-front thick- 
ness is expected to grow over the lifetime of a quasar given 
the discrepancy in evolution timescales between the largest 
and smallest scales. At small radii, photo-ionization equilib- 
rium is reached quickly since ionizing photons are abundant, 
whereas geometrical dilution and attenuation of the initial ra- 
diation field slow ionization evolution considerably on large 
scales, effectively 'stretching out' the I-fronts of hydrogen 
and helium with time. A monochromatic representation of the 
quasar SED leads to a reduction in this effect, but also leads 
to severe errors in the overall ionization structure (see Figure 
3 of KH08). These errors are of the same order of magnitude 
as those resulting from the neglect of physical effects, such as 
ionization via helium recombination photons (KH08, Figure 
6), or ionization from secondary electrons (KH08, Figure 7). 
These effects are likely important in studies of radiative feed- 
back from stars and active galactic nuclei (AGNs), and most 
certainly in efforts to simulate cosmological reionization. An 
effort must be made to ensure that the SEDs used in numerical 
simulations accurately reflect the properties of their continu- 
ous analogs, especially if it is spectrum-dependent effects in 
which we are most interested. 

We will focus on the following questions in this paper. How 
significant are the errors in the temperature and ionization 
state of a medium that arise solely due to the discretization of 
SEDs? How many frequencies are required to minimize such 
errors, where must they be positioned in frequency-space, and 
how should their relative luminosities be apportioned? For 
what numerical methods is it possible to represent sources 
with continuous SEDs, or are there perhaps advantages in dis- 
cretizing SEDs, even when it is not required by the algorithm 
of choice? Answers to these questions may lead to revised 



interpretations of previous studies which used discrete radia- 
tion fields, but more importantly, will reduce the guesswork 
involved in discretizing SEDs, and promote frequency resolu- 
tion to the same status as spatial, temporal, and mass resolu- 
tion, which are more easily selected on a problem-by-problem 
basis. 

In Section [2] we will introduce the one-dimensional ra- 
diative transfer framework used to obtain the solutions pre- 
sented in later sections. In Section [3] we quantitatively as- 
sess the accuracy with which multi-frequency calculations 
reproduce the ionization and heating profiles of continuous 
SEDs. Section @] is devoted to introducing a technique for 
optimally selecting discrete SED templates, and Section [5] 
will present the results obtained with this method, including 
applications to one-dimensional and fully three-dimensional 
radiation-hydrodynamic calculations. Discussion and conclu- 
sions can be found in Sections [6] and [7] respectively. Valida- 
tion of the radiative transfer code used for this work and fur- 
ther details regarding the optimization algorithm can be found 
in the Appendix. 



2. RADIATIVE TRANSFER FRAMEWORK 

One dimensional radiative transfer calculations around 
point s ources have been used to mo del cosmological reion- 
ization (Fukugi ta & Kawasakill 19941) . the thickness of quasar 
ionization fronts (KH08), the time-evolution of ioniza- 
tion and heating around first stars, galaxies, and q uasars 
(iThomas & Zaroublll2008t IVen katesan & Bensonll201 ll) . and 
their associated observable signatures. Given that our focus is 
on frequency resolution, it would be unnecessary to perform 
calculations in a more complex setting than this, with addi- 
tional unrelated physics. As a result, our one-dimensional 
methods strongly resemble those used by previous authors, 
though for completeness, we will reiterate the aspects of these 
methods most pertinent to the problem at hand. 

In general, the chemical and thermal evolution of gas sur- 
rounding a radiation source is governed by a set of differential 
equations describing the number densities of all ions and the 
temperature of the gas. Assuming a medium consisting of hy- 
drogen and helium only, we first solve for the abundances of 
each ion via 



dn Hu 

dt 
dn He n 

dt 



driHe 



(r H I + YH 1 + Ph I»e)«H I - «H n«e«H II (1) 
(rHe I + YHe I + Pile I«e)«He I + OCHe III«e«He III 



dt 



(pHe II + CtHe II + ^HeIl)"e"He II 
: (rHe II + YHe II + pHe II«e)«He II ~ OCHe m«e«He 



(2) 



(3) 



Each of these equations represents the balance between ion- 
izations of species H I, He I, and He II, and recombinations 
of H II, He II, and He III. Associating the index i with ab- 
sorbing species, i =H I, He I, He II, and the index 1 with ions, 
i' =H II, He II, He III, we define r, as the photo-ionization 
rate coefficient, Y as the secondary ionization rate coefficient, 
a,/ as the case-B (dielectric) recombination rate coeffi- 
cients, p,- as the collisional ionization rate coefficients, and 
«c = «H ii + «He n + 2«He hi as the number density of electrons. 
At each time step, we also solve for the temperature evolu- 
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tion, dTk/dt, which is given by 

lit (^T*) =^%nM-%^-%Wi> 

- ¥i»e«i - C0 He n w e«He „ (4) 

i 

where J-(i is the photo-electric heating rate coefficient (due 
to electrons previously bound to species i), <x>H e n is the di- 
electric recombination cooling coefficient, and t[f, and 
\|// are the collisional ionization, recombination, and colli- 
sional excitation cooling coefficients, respectively. The con- 
stants in Equation (|4]i are the total number density of baryons, 
«tot = "H + «He + «e, the mean molecular weight, /j, Boltz- 
mann's constant, k%, and the fraction of secondary electron 
energy dep osited as heat, f heat . We use the formulae in Ap- 
pendix B of Fukugita & Kawasaki ( 19941) to compute the val- 
ues Of a,, P,-, Ci. T li'. Vi. and ro He II- 

The most critical aspect of propagating the radiation field in 
our one-dimensional simulations is computing the ionization 
(r,, y) and heating (9-Q) rate coefficients accurately. In order 
to directly relate our results to fully three-dimensional radia- 
tive transfer calculations, we have chosen to adopt a photon- 
conserving (PC) algorithm nearly identical to those employed 
by several widely used c odes, like C 2 Rav (IMellema et alj 
l2006t |Friedrich etaDHoH), and Enzo dWise & Abelll201 11) . 
Our code is able to compute T,, y, and ^ in a non-photon- 
conserving (NPC) fashion as well, to ena ble comparison with 
previo us one-dimensional work such as Thomas & Zaroubi 
(2008). The two formalisms are equivalent in the limit of 
very optically thin cells, a condition that can be met easily 
in one-dimensional calculations but is rarely computationally 
feasible in three dimensions. For NPC methods, if the opti- 
cal depth of an individual cell is substantial, the number of 
ionizations in that cell will not equal the number of photons 
absorbed for that cell, i.e., photon numb er will not be con - 
served. This problem was remedied by lAbel et all (1999), 
who inferred the number of photo-ionizations of species i in a 
cell from the radiation incident upon it and its optical depth, 

At,.v = rijOj.yAr. (5) 

It is most straightforward to imagine our one-dimensional 
grid as a collection of concentric spherical shells, each having 
thickness Ar and volume V S h(r) = 4n[(r + Ar) 3 — r 3 ]/3, where 
r is the distance between the origin and the inner interface of 
each shell. The ionization and heating rates can then be re- 
lated to the number of absorptions in any given shell (thus 
preserving photon number), as 

H=A,J\y-y i )I v e-^ (\-e-^) ^, (8) 

where we have defined the normalization constant A, = 
Lboi/n/VshfV), and denote the ionization threshold energy for 
species i as /zv,. / v represents the SED of radiation sources, 
and satisfies Ll v dv = 1, such that LboiA' = Lv- 

Equation (01 represents ionizations of species i due to fast 
secondary electrons from photoionizations of species j, which 



has number density nj, and ionization threshold energy, hVj. 
jf n is the fraction of photo-electron energy deposited as ion- 
izations of species i. In the remaining sections we only in- 
clude the effects of secondary electrons when considering X- 
ray sources, which emit photons in the range 10 2 eV < E < 
10 4 eV. In this regim e, the values of f heat and f}° n computed 
via the formulae of Shul l & van Steenb erg (1985) are suffi- 
ciently accurate, but for radiation at lower energies where 
yheat an( j y-ion h ave a stronger ener gy dependence, the fit- 
ting formulae of | Ricotti et all (120021) or the lookup tables of 
Furlanetto & Stoeverl (120101) would be more appropriate. The 
total secondary ionization rate for a given species, y,, is the 
sum of ionizations due to the secondary electrons from all 
species, y = ff a y, 7 n } j n t . 

The optical depth, T v = T V (V), in the above equations is the 
total optical depth at frequency v due to all absorbing species, 
i.e., 

T v(r) = V f Giyni{r')dr 

= £a iiV AT / (r) (9) 

i 

where Nj is the column density of species i at distance r from 
the source. We calculat e the bound-free absorption cross- 
sections using the fits of [Verner et al.l (fl996) throughout. 

The values of T,, y , and are completely predetermined 
for a given radiation source, and as a result, can be tabulated as 
a function of column density to avoid evaluating the integrals 
in these ex pressions numerically 'on-the-fly' as a simulatio n 
runs (e.g.. IMellema et al.1 120061 iThomas & Zaroubil 120081) . 
Isolating the frequency-dependent components of Equations 
we can define the integrals 

r°° dv 
<t>,(x v )= / V~ Xv — (10) 
Jv, nv 

Vj(l\r) = (°°Ive-^dv, (11) 
J\i 

allowing us to re-express the rate coefficients as 
r,=A,[4>,(^v)-4>,(<v)] (12) 

Yy = ^ {^N) - ^K-, v ) - kvj [4> ; -(Tv) - <&;(t} >v )] } 

(13) 

H = A, {Wi{z v ) -W/«- v ) -hVi [*i(Tv)-*i(i/ |V )] } . (14) 

where = T v + Ax,-, v - Later references to "continuous 
SEDs" signify use of this technique, where the integral val- 
ues <!>,■ and ^P, are computed over a column density interval of 
interest a priori using a Gaussian quadrature technique, rather 
than on-the-fly via discrete summation. 

Tabulating Equations ([Tot and (JTTJ grants a significant 
speed-up computationally, but also forms the basis of our fre- 
quency resolution optimization strategy (Section H). Note, 
however, that in general the dimensionality of these lookup 
tables is equal to the number of absorbing species (through 
Ax, )V ), so the tables for simulations including hydrogen 
only are one dimensional, while those including hydrogen 
and helium are three dimensional. If we chose to ad opt 
the secondary electron trea tment of iRicotti et al.l (12002b or 
Furlanetto & St oeveri (120101) . our lookup tables would inherit 
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an additional dimension, as the secondary ionization and heat- 
ing factors fl° n and / heat would depend both on photon energy 
and the hydrogen ionized fraction, xh u. 

Equations (TT2T> — <fT4-b are completely general for PC algo- 
rithms, whether the source SEDs are discrete or continuous 
— the only difference being for discrete SEDs, the integrals 
in Equations ( fTOb and ( fTTT ) become sums over the number 
of discrete emission frequencies, n v . In practice, computing 
r,, Ji, and % is more straightforward for sources with dis- 
crete SEDs, as we can simply count the number of ionizations 
caused by photons at each individual frequency, and convert 
this into the amount of excess electron kinetic energy avail- 
able for further heating and ionization. When testing the ac- 
curacy of discrete solutions in later sections we employ this 
method, where radiation is emitted at n v frequencies, with 
each frequency v„ carrying a fraction /„ of the source's bolo- 
metric luminosity. The photoionization and heating coeffi- 
cients can then be expressed as 

r,,„ = ^V x Mi-^ Ax '' v '') d5) 

hv„ 

7y> = r Wv„-V;)/v/ (16) 

^ I „ = r /i v>(v B -V i ). (17) 

The total rate coefficients can be found by summing each 
of these expressions over all frequencies, n = 1,2,3, ... ,n v . 
These equations are identical to Equations (fT2l-([T4li for the 
discrete SED case, but are perhaps more intuitive. 

For simplicity, our current treatment neglects a few physi- 
cal processes that are cosmological in origin, or simply do not 
rely on the radiation field directly. These include cooling via 
free-free emission and hydrogen and helium ionization due 
to helium recombination photons (which depend on the gas 
kinetic temperature and electron density), and cosmological 
effects such as Hubble cooling, Compton cooling off cosmic 
microwave background (CMB) photons, and photo-ionization 
by Wien-tail CMB photons (which depend on kinetic temper- 
ature, redshift, and the Hubble parameter). 

Two additional approximations are implicit in the remain- 
der of this paper. They are (1) the infinite speed-of-light 
approximation and (2) the on-the-spot approximation (we 
use the case-B recombination coefficients in Equations (Q]i- 
(O). The former approximation could be dubious for very 
bright sources in low-density media, while the latter is gen- 
erally not a good as sumption, as discussed at length in 
Cantalupo & Po rcianil (1201 ll) . As a result, the absolute accu- 
racy of our solutions is not guaranteed in regimes where care - 
ful treatment of the speed of light and recombination photons 
is necessary, but this is acceptable since we only care about 
the relative differences among our solutions. The optimized 
SEDs of Section [5] will apply equally well to simulations in- 
cluding more ionization and/or heating/cooling processes, so 
long as they do not depend directly on the radiation field (e.g., 
ionization of H I and He I by helium recombination photons; 
iFriedrich et aDl2012h . 

3. ASSESSING THE CONSEQUENCES OF DISCRETE RADIATION 

FIELDS 

To quantify the differences between the ionization and tem- 
perature profiles around sources with continuous and discrete 
SEDs, we will simulate two test problems. First, the standard 
case of a 10 5 K blackbody in a hydrogen-only medium, and 
second, a power-law X-ray source in a medium consisting of 
both hydrogen and helium. 



3.1. 10 5 K Blackbody 

The 10 5 K blackbody problem has been studied extensively 
(e.g., T est Problem 2 in the Radiative Transfer Comparison 
Project; Ili ev et al.l200 6, hereafter RT06) due to its simplicity, 
and perhaps also because the surface temperatures of PopIII 
stars are expected to be ~ 10 s K (Schaerer 2002). We adopt 
nearly the identical setup as in RT06, i.e., a uniform hydrogen- 
only medium with number density «h = 10~ 3 cm~ 3 , ini- 
tial ionized fraction ihh = 1.2 x 10~ 3 , initial temperature 
T = 10 2 K, and a 10 5 K blackbody with an ionizing pho- 
ton luminosity of Q = 5 x 10 48 s _1 . The only difference 
between our simulations and RT06 is that we use a domain 
Lbox =10 kpc in size, rather than L\, ox = 6.6 kpc, to allow for 
a comparison of discrete and continuous solutions at slightly 
larger radii. We evolve the simulations for 500 Myr on a grid 
of 200 linearly spaced cells between 0. 1 < r/kpc < 10, ignor- 
ing the details of secondary ionization (i.e., all photo-electron 
energy is deposited as heat). 

In Figure Q] we compare the ionization and temperature 
profiles around two 10 5 K 'blackbody' sources of constant 
ionizing photon luminosity Q = 5 x 10 48 i _1 — one a true 
blackbody emitter with a continuous SED spanning the range 
13.6-100 eV (black lines), and the other with a monochro- 
matic SED at hvi = 29.6 eV, the average energy of ioniz- 
ing photons for this source (red lines). We can see the same 
qualitative results that have been pointed out by previous au- 
thors, namely, that monochromatic sources of radiation fail to 
ionize (top panels) and heat (lower panels) gas at large radii 
as significantly as continuous sources, since all photons are 
absorbed near a single characteristic column density, repre- 
senting the point where x V] « 1, i.e., ./V c har ~ Cy/- Th e re l a ~ 
tive error in the position of the ionization front, Arip, where 
rjF = r{xn i = x H „ = 0.5), is 8% after 10 Myr, 10% after 100 
Myr, and 1 1 % after 500 Myr. In the optically thin regime, the 
monochromatic spectrum overestimates ionization by factors 
of two to three on average and up to an order of magnitude at 
all times, though the latter effect is primarily because the neu- 
tral fraction is a steeply declining function with decreasing 
radius, and the I-fronts of the two solutions are offset. Out- 
side the I-front, the situation is more interesting as the gas is 
mostly neutral. After 100 Myr of evolution, the ionized frac- 
tion outside the I-front is underestimated by a factor of two on 
average, and by as much as a factor of six. 

The temperature evolution, shown in the bottom panels of 
Figure Q] is significantly more troubling. The monochro- 
matic source captures the temperature well within the ioniza- 
tion front where the gas is in photoionization equilibrium, but 
quickly diverges from the continuous solution outside. Like 
the ionization profiles, discrepancies grow with time. After 10 
Myr of evolution, the monochromatic source underestimates 
the temperature at large radii by a factor of two on average, 
and by a factor of seven at the point of greatest discrepancy. 
After 100 (500) Myr, the discrete solution underestimates the 
temperature by up to a factor of 17 (41). 

If considering the heating and ionization around a single 
PopIII star, the errors induced by monochromatic treatments 
may not be cause for concern upon first inspection since 
PopIII stars are expected to live only a few Myr, and we can 
see that errors are less significant at early times. However, the 
intergalactic medium (IGM) is subject to the ionization and 
heating caused by all sources, whose cumulative impact will 
be substantial even though the ionization and heating caused 
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t = 100 Myr 




Figure 1. Comparison of ionization (top) and temperature (bottom) profiles around a 10 5 K blackbody source after 10 Myr (left) and 100 Myr (right) using 
continuous (black) and monochromatic (red) SEDs. Solid lines in the top panels correspond to the neutral fraction (ah [), while dashed lines correspond to the 
ionized fraction (xu n). We apply these line color and line style conventions for all radial profiles presented in this paper. 



by individual sources may be very small. Globally, then, the 
IGM is insensitive to individual stellar lifetimes, and instead 
evolves as it would if ionizing photons originated from a sin- 
gle, very luminous, very long lived object. 

This manner of thinking has already materialized in the 
realm of large volume cosmological simulations, where 'star 
particles' are generally as luminous as one or more star clus- 
ters, and 'galaxy particles' behave in a way that is consistent 
with the integrated properties of an entire galactic stellar pop- 
ulation (and perhaps active nucleus). Such approximations are 
necessary with limited spatial resolution, but more than ade- 
quate for studies of the IGM. Over time though, errors in gas 
properties due to poor frequency resolution will accrue, as it is 
the combined properties of all radiation sources which affect 
IGM properties, however short-lived each individual source 
may be. 

3.2. Power-Law X-Ray Source 

To address the effects of discrete SEDs in environments 
where multiple chemical species are important and large at- 
tenuating columns are possible, we now turn our attention 
to a power-law X-ray source embedded in a 1 Mpc domain 
consisting of hydrogen and helium, with a primordial helium 
abundance (by mass) of Y — 0.2477. 



Our selection of parameters for this problem is motivated 
by studies of high-redshift quasar s, and particularly their role 
in the epoch of reionization (e.g., Venkatesan et al. 200j]). X- 
rays have long mean free paths, and as a result are capable of 
ionizing and heating gas on very large (^Mpc) scales. Large- 
scale heating is responsible for driving the high-redshift all- 
sky 21 cm signal toward emission, and inducing fluctuations 
in 21 cm power spectra on large angular scales (for a review of 
21 cm cosmology, see Furlan etto et al.l (12006b ). An early X- 
ray background may also be important in interpret ing the opti- 
cal depth to electron scattering o f the CMB (e.g jRicotti et al.l 
l200HIShull & Venkatesanll2008l) . 

While supernovae and/or X-ray binaries could be impor- 
tant sources of hard photons in the early universe, we as- 
sume the source of X-rays is persistent — an accreting SMBH 
with mass M. — 1O 6 M and radiative efficiency of e. = 10%, 
which leads to a bolometric luminosity of L\, ] = e.ilgdd — 
1.26 x 10 43 erg s -1 . Here, Z^d = 4nGM.m p c/aT is the Ed- 
dington luminosity, where m p is the proton mass and Ot the 
Thomson cross-section. The mass (and thus luminosity) of 
the SMBH is allowed to grow as it accretes, 



M,(/)=M,(0)exp 



1-e. 
e. 



t 

fedd 



(18) 
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where f e( jd = 0.45 Gyr is the e-folding timescale for SMBH 
growth (an Eddington, or Salpeter time). The SED is taken to 
be a power law of the form 



IkeV J 



1-oc 



(19) 



where a is the spectral index. We adopt a = 1.5, over the en- 
ergy range 10 2 -10 4 eV. The surrounding medium has a con- 
stant mass density of p = 5.4 x 10~ 28 g cm~ 3 (cosmic mean at 
redshift z= 10), initial ionized fractions xh n = XHe n = 10~ 4 , 
*He m = 0, and initial temperature Tq = 10 2 K. The domain for 
this problem is divided into 400 cells linearly spaced between 
0.01 < r/Mpc < 1, and is evolved for £.f e dd = 45 Myr. 




Figure 2. Comparison of hydrogen (top) and helium (bottom) ionization pro- 
files around an a = 1.5 power-law X-ray source after 45 Myr using continu- 
ous (black) and monochromatic (red) SEDs. 



In Figure [2] we compare the hydrogen and helium ioniza- 
tion profiles for two X-ray sources having the same bolo- 
metric luminosity. One, a continuous power-law source as 
described above, and the other a monochromatic source of 
0.5 keV photons (a fiducial monochromatic emission energy). 
The monochromatic source underestimates the radii of both 
the hydrogen and helium ionization fronts by a factor of ~ 2.3, 
and overestimates the hydrogen neutral fraction on average by 




Figure 3. Comparison of temperature profiles around an a = 1.5 power-law 
X-ray source after 45 Myr using continuous (black) and monochromatic (red) 
SEDs. 



a factor of three, and at most by a factor of 20 within the hy- 
drogen I-front. The same general picture applies to helium, 
where errors in the neutral helium fraction are enormous since 
the He I-He II I-front is very sharp (as it was for hydrogen in 
the previous section), and xn e n and xn e m are in error by fac- 
tors of 2-20 depending on radius. 

Errors in the temperature profile are less extreme, as shown 
in Figure [3] On small scales, the monochromatic source 
captures the temperature quite well, but at large radii, the 
monochromatic source overestimates temperatures by a fac- 
tor of two on average. 

The disparity in the magnitude of ionization and tempera- 
ture errors is a reflection of the strong frequency dependence 
of the bound-free absorption coefficients. Photo-ionization of 
hydrogen or helium by 0.5 keV photons is rare, but when it 
does occur, at least ~ 90% of the original photon energy is left 
to be deposited mostly as heat, unless the free electron density 
is very low. Because the ionization of hydrogen and helium 
by the monochromatic source is very inaccurate, errors in the 
free electron density will substantially alter the amount of sec- 
ondary electron energy deposited as heat, rather than further 
ionization. 

The consequences of miscalculating ionization and heating 
could affect efforts to model and interpret current and future 
21 cm measurements, since the primary 21 cm observable, the 
differential brightness temperature (87],), depends on the hy- 
drogen neutral fraction, UV radiation field, electron density , 
and the gas kinetic temperature (Tk) dFurlanetto et al.l 2006). 
Neglecting the presence of a Lya background, the scaling 



5T b o,T^(l+8)(l+z)- 1/2 x\T ne 



, »e > «H I 
, «e < «H I 



(20) 



holds approximately in regimes where 7cmb "C 7k < 10 4 K. 

In the immediate vicinity of radiation sources where gas is 
entirely ionized, 87], — > due to the leading xh i term, but at 
large radii where the ionizing flux is weaker, the 87], signa- 
tures of stars and quasars could vary significantly solely due 
to miscalculations of xu u n e , and Tk- The above scalings have 
especially strong consequences for gas within a few Mpc of 
strong X-ray sources, where hydrogen is weakly ionized, tem- 
peratures are of order 10 2 -10 3 K, and the free electron density 
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is enhanced due to efficient ionization of helium by the hard 
radiation field. In the earliest stages of reionization where 
Tk < Tcmb (z) and the Lya background is important, errors in 
jch i, «e, and Tk will lead to errors in Sli as well, though in a 
less straightforward way, since the spin temperature, T$, must 
be computed carefully. 



4. OPTIMIZATION STRATEGY 

To avoid errors of the sort described in the previous section, 
we have developed a technique for optimally constructing dis- 
crete SEDs that preserves the ionization and heating proper- 
ties of their continuous counterparts. Although ray-tracing 
algorithms are capable of tabulating the relevant ionization 
and heating quantities (Equations dlOb and (fTTT i). few codes 
have taken advantage of this, and have instead cast monochro- 
matic r ays (e.g., state of the art reionization simulations with 
wv = 5:lTrac et al.ll2008l) . Monte Carlo codes (e.g., CRASH; 
iMaselli et al.l 120031) have been used to simulate reionization 
with w v > 20 multi-frequency photon packets (Ciardi et al. 
120121) . though such a large number of frequencies may be 
computationally debilitating for some algorithms, or unnec- 
essary depending on the problem of interest. 

Even when the algorithm of choice is compatible with prop- 
agating continuous radiation fields via tabulation of Equations 
( flCTb and ( fTTI ). it may not be computationally advantageous. 
The overhead alone can in fact be substantial, particularly in 
the case of source-dependent SEDs — for example, the SED 
of a stellar population as a function of age, or BH accretion 
spectra that vary with mass or luminosity. Such situations 
would require a separate lookup table for Equations (flOl and 
( fTTT i at each age/mass/luminosity of interest for a given radia- 
tion source. In addition, there are algorithms for which prop- 
agating continuous radiation fields in large volumes become 
completely intractable, yet large volumes are a necessity for 
the science questions of interest (e.g., reionization). For more 
discussion on these issues, see Section|6] 

As introduced in Section|2] our optimization strategy relies 
on the fact that the SED of a radiation source appears only in 
the quantities <t>, and *Pf (see Equations ( fTob and (fUJi). If we 
can construct a discrete SED that reproduces the values of <!>, 
and *P/ to a high degree of accuracy over a column density 
interval of interest, then the discrete radiation field is indistin- 
guishable from its continuous counterpart, and we have suc- 
cessfully preserved the true radiative properties of the source. 

For sources with discrete SEDs, Equations (flOl and (fill 
become 

n v J 

#tM = I, Pi) 

" v 

?!W = EV^, (22) 

n=l 

where we have used primes to indicate that these quantities 
are computed by direct summation over n = 1,2, . . . ,n v fre- 
quencies, rather than by a continuous integral. 

Ensuring that <t>, = <£>| and *P; = *F| is a minimization prob- 
lem of dimensionality 2n v , since each additional frequency 
bin lends two degrees of freedom — its frequency (v„), and 
the fraction of the bolometric luminosity assigned to that fre- 
quency (/„). Our goal is to minimize the difference between 



continuous and discrete solutions, i.e., 

4> ; - = 

¥,-¥• = 0. (23) 

These functions span several orders of magnitude over a broad 
range in column density, making it more practical to seek so- 
lutions to 



log(|j)=0 (24) 

which place equal emphasis on all column densities. Pre- 
serving the high column density behavior of <£>,■ and is 
especially important for very luminous sources and/or en- 
vironments with dense clumps in the immediate vicinity of 
the source, since the actual photoionization and heating rates 
are a combination of <£>,-, VP,-, and the normalization factor 
Ai oc L bol /r 2 . 

For a given n v and source SED, we solve Equation 
(|24"1 > using the optimization technique Simulated Annealing 
dKirkpatrick etalj fl983: Cerny 1985), which tr averses our 
2n v dimensional parameter space in search of the frequency- 
normalization pairs (v„,/„) that best reproduce the values of 
<!>; and ^P, . We leave a more detailed description of the algo- 
rithm and our implementation of it to the Appendix. 

5. RESULTS 
5.1. Optimal Discrete SEDs 

We have obtained optimal SEDs for a 10 5 K blackbody 
emitting in the range 13.6-100 eV, and ana= 1.5 power- 
law X-ray source with emission spanning the interval 10 2 - 
10 4 eV. In each case, we set the upper column density limit 
for our optimization to be the column density of a fully neu- 
tral medium, i.e., f = nH^box and A^ff = «He£box, where 
we use Lbox to denote the size of the domain, as in RT06. 
For the 10 5 K blackbody simulations, this works out to be 
N™ x = 3.1 x 10 19 cm~ 2 , and for the power-law X-ray simula- 
tions, N™ x ~ x 10 22 cirr 2 and A^ e a * - x 10 21 cm" 2 For cos- 
mological simulations with periodic boundary conditions, the 
upper column density limits would need to be chosen based on 
a maximum length scale of interest, or for radiative feedback 
focused simulations, by the column density of the densest ob- 
jects of interest (damped Lya systems, for example). Such 
choices are already made in ray-tracing calculations to limit 
computational expense. Generally, rays are terminated once 
the emission has been attenuated by a large factor. 

The only situation in which we do not evaluate the full cost 
function is ra v = 1> where we instead optimize for the opti- 
cally thin regime alone (i.e., only the first term of Equation 
IA2I ). where <t>, and ^ are ~ constant with column density. In 
this case, the optimal solutions are simply those that preserve 
the bolometric luminosity of the source and the total number 
of ionizing photons, and can be verified analytically (Equa- 
tions ( [Tol l and (fTTT i). For the case of a hydrogen and helium 
medium, we have found that neglecting He II opacities miti- 
gates the computational cost of the computation while result- 
ing in no appreciable changes in our optimal SEDs and thus 
negligible changes in <t>' and *F'. The main results are summa- 
rized in Figures|6]and|7]and Tables[T]and[2] all results derived 
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Table 1 

Optimal SEDs for 10 5 K Blackbody Sources 





n=\ 


n = 2 


n = 3 n = 4 


1 


(29.61.0.89) 






2 


(27.93,0.68) 


(62.04,0.21) 




3 


(20.58,0.39) 


(40.75,0.39) 


(69.23,0.11) 


4 


(17.98,0.23) 


(31.15,0.36) 


(49.09,0.24) (76.98,0.06) 



Note. — Each entry is the (h\ n ,I n ) pair for bin n. Energies 
are in units of eV, and normalizations are expressed as fraction 
of the bolometric luminosity. 



Table 2 

Optimal SEDs for a = 1.5 Power-Law X-ray Sources 



»v 


n = \ 


n = 2 


n = 3 n = 4 


1 


(999.98.1.00) 






2 


(255.87,0.17) 


(2553.6,0.83) 




3 


(171.93,0.08) 


(518.22,0.14) 


(3098.5,0.78) 


4 


(146.11,0.05) 


(307.30,0.07) 


(704.56,0.14) (3564.2,0.73) 



Note. 

source. 



Same as Table Q]but for an a = 1.5 power-law X-ray 



from K = 2 x 10 4 and K = 10 4 Monte-Carlo trials, for the 10 5 
K blackbody and a = 1 .5 power-law source, respectively. 

From Tables Q] and [2] it is clear that the optimal emission 
frequencies for both sources are not evenly spaced above the 
hydrogen or helium ionization thresholds, either in linear or 
log-space. In each case, the addition of a new frequency bin 
leads to a decrease in both the emission frequency and nor- 
malization of all other bins. This signifies (1) the efficacy with 
which high energy photons photoionize and photoheat gas at 
large column densities (a regime inaccessible to lower energy 
photons which become optically thick at small columns), and 
(2) the increase in excess electron kinetic energy available for 
further ionization and heating with increasing photon energy. 
The former effect is most important for the blackbody source, 
which we can see in Figure [4] Not surprisingly, it is the low- 
est energy photons {hv\ = 17.98 eV) in the n v = 4 spectrum 
that are responsible for the ionization (through <t>) in the op- 
tically thin regime, while successively higher frequency bins 
become the primary agents of ionization as we move to higher 
column densit ies. The same trend does not hold completely 
in Figure [4(b)] as in this case it is the second and third energy 
bins that provide the bulk of the heating (through *P) at low 
column densities. 

For the X-ray source, the second effect dominates, as the 
optical depth at any column density is small for most photons 
considered (10 2 < hv < 10 4 eV) over the entire domain. As 
shown in Figure[5] the photons responsible for the majority of 
the heating (through *P) over all column densities are those in 
the highest energy bin, the same photons which are the least 
effective at ionization. The trends and errors of Figure |5]are 
the same for <t>, and *P as a function of helium column density. 

In Figures [6] and [7] we show the probability distribution 
functions (PDFs) for the position and normalization of the op- 
timal SED frequency bins obtained (drawn from Tables [T]and 
ID). Solutions are less tightly constrained as n v is increased, as 
evidenced by a broadening in the distributions of frequency 
and normalization for each bin. This behavior is expected, 
given that each new bin contributes to the magnitude of <t> 
and *P in some region of column density space previously oc- 



K-H-H-H 





Figure 4. Top Panels: Comparison of 3>h i and Q>' a , (a) and "Ph i and *Py , 
(b) as a function of H I column density for a 10 5 K blackbody, showing the 
numerically computed continuous integral (solid black), best-fit composite 
four-bin discrete sum (blue crosses), and the contribution from each individ- 
ual discrete frequency bin (dashed blue). Annotations represent the (/iV„,/„) 
pah's for each frequency group, drawn from Table [T] Bottom Panels: Percent 
error between discrete and continuous solutions. The solid blue line is the er- 
ror for the four-bin optimal solution, while the errors induced by three-, two-, 
and one-bin solutions are shown in magenta, green, and red, respectively. 

cupied by one or more other frequencies. 

Holding /„ constant, a decrease in v„ will cause a negative 
vertical shift in the contribution of bin n to the magnitude of 
<t>, for example, but will simultaneously add power at larger 
column densities, since the turnover point for bin n occurs 
at A'char ~ cr^ 1 , and a v „ ~ V~ 3 . To avoid an increase in /, 
the power lost at small column densities has to be compen- 
sated for, either by a decrease in V„_i, or an increase in /„_i, 
where n — 1 denotes the bin with frequency v„_i < V„. As a 
result, there are degeneracies between all bins, and the mag- 
nitude of the degeneracy is greatest for bins positioned closest 
in frequency-space. In order to tighten the PDFs for each opti- 
mal frequency bin, one or more terms would need to be added 
to /, in order to assign preference to one set of bins over an- 
other. For our purposes, any SED that minimizes / is just as 
good as any other, but additional terms in the cost function 
are certainly justifiable in the case of a ray-tracing calcula- 
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Figure 5. Same as Figuref4]but for an a = 1.5 power-law X-ray source. 

tion, where higher emission frequencies increase the compu- 
tational cost of a calculation since their mean free paths are 
long. Adding a term to / that scales with v„ would encourage 
optimal SEDs with the smallest emission frequencies possi- 
ble, for example. 

Optimization for n v > 4 is certainly possible, though un- 
necessary in our case. At a given frequency, the transition 
from optically thin x = to optically thick (x > 1) in the 
functions <t> and SP occurs over an order of magnitude in 
column density (by definition, see Equation (0). For both 
SEDs we have investigated, the column density regime of in- 
terest spans fewer than four orders of magnitude, motivating 
our choice of 1 < « v < 4. We have performed optimizations 
with n v > 4, but the addition of each additional bin when 
n v > logjoC^Vmax/^Vmin) reduces the error between <t> and <!>', 
and *P and *F' much less significantly than additional bins 

when « v < log 10 (AW x /Af mi n). For a § iven "v. increasing N ma x 
will simply increase max|<S> — <£>'| and max^ — \P'|. 

5.2. Confirmation with One-dimensional Calculations 

To verify the solutions of the previous section, we ran sim- 
ulations identical to those of Section [3] but with our optimal 
discrete SEDs. We compute T,, y,, and via Equations 
(TT3T) — (TT7l> "on-the-fly," rather than generating lookup tables 
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Figure 6. Emission energy (a) and normalization (b) probability distribution 
functions (PDFs) of optimized discrete 10 5 K blackbody spectrum using n v = 
1,2,3,4 (from bottom to top). In each panel, the gray histogram denotes the 
initial guesses for all Monte-Carlo trials, and the black, blue, red, and green 
histograms show the end point for the first, second, third, and fourth bins, 
respectively (ordered by increasing emission frequency). 

of <!>, and M*,-. As expected, accurate preservation of the quan- 
tities 4>, and ^P, over the column density ranges of interest 
renders ionization and temperature profiles around sources 
of discrete radiation indistinguishable from their continuous 
counterparts. 

In Figure [8] we compare ionization and heating around a 
10 5 K blackbody after 100 Myr of evolution as in Section [3j 
showing the solution obtained with our optimal monochro- 
matic (red) and four-bin (blue) SEDs. The continuous and 
four-bin solutions are indistinguishable. 

In Figure [9] we perform the same analysis for the a = 1 .5 
power-law simulations. Our optimal four-bin SED reproduces 
the hydrogen and helium ionization profiles (and thus electron 
density) and temperature of a continuous SED to high preci- 
sion. The most noticeable errors are in the hydrogen neutral 
fraction within the hydrogen ionization front, where errors be- 
tween four-bin and continuous solutions are still only ~ 1%. 
Errors in xn e m are negligible, justifying our neglect of A^e n 
in the optimization process. 

It should be noted that our optimal monochromatic SED 
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Figure 7. Same as Figure[6]but for an a = 1.5 power-law X-ray source. 

for the X-ray source performs even more poorly than the 
fiducial 0.5 keV SED. This signifies a general problem with 
monochromatic emission for any spectrum with a hard com- 
ponent. Whereas the monochromatic optimization (t v = 0) 
works quite well in the 10 5 K blackbody case since hydro- 
gen absorbs UV photons readily, X-rays are not so readily 
absorbed by hydrogen and/or helium. As a result, the char- 
acteristic column density where most 1 keV photons are ab- 
sorbed lies outside of our domain, leading to severe under- 
ionization (of all species) and under-heating. The reason the 
0.5 keV SED works better is because its characteristic absorp- 
tion column is smaller, lying within our domain. We have 
experimented with relaxing the optically thin requirement for 
monochromatic optimization, and find that it is equally diffi- 
cult to preserve ionization and heating profiles with emission 
at a single frequency. 

5.3. Three-dimensional Radiation-hydrodynamic Simulations 

with Enzo 

To study the impact of spectral discretization in a more 
complex setting, we ran RT06 test problem 2 with hydro- 
dynamics, as well as two fully three-dimensional cosmolog- 
ical radiation-hy drodyna mic simulati o ns sim ilar to those of 
lAbel et al.l (12007b and lAlvarez et all <|2009), both with the 




g 10 

% 10 3 



10" 
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Figure 8. Comparison of ionization (top) and temperature (bottom) profiles 
around a 10 s K blackbody source after 100 Myr showing the solutions ob- 
tained using continuous (black), monochromatic (red), and optimal four-bin 
discrete (blue circles/squares) SEDs. 



Enzo code dBrvan & Normanll997tlO'Shea et al ll2004lf T l. AU 
analysis was performed with yt dTurk et al.ll201 ll) . 

The results of the RT06 radiation-hydrodynamic test prob- 
lem are shown in Figure [10] where we compare the solutions 
obtain ed using the four-bin SED employed by IWise & Abell 
(201 1) in addition to our own (TableHJ. The solutions are in- 
distinguishable, which is expected given the relatively small 
range of column density explored in this problem. 

The cosmological simulations follow the formation of a 
100M© PopIII star, its brief 2.7 Myr lifetime in which it 
emits 1.2 x 10 50 ionizing photons per second, and the X- 
ray emission resulting from accretion onto a remnant BH as- 
sumed to form via direct collapse after stellar death (as in 
lAlvarez et al.l (120091) ). The accretion rate, and thus luminos- 
ity assuming e. = 10%, is the Bondi-Hoyle accretion rate of 
the cell in which the BH resides. The simulation volume is 
0.25 Mpc h on a side, with 128 3 particles and cells on the 
root grid. A single nested grid occupies the inner 1/8 of the 
volume at twice the root grid resolution, where eight addi- 
tional levels of adaptive-mesh refinement are allowed, yield- 
ing a peak spatial resolution of 0.23pc h . 

1 Revision f 4a8b5f 5e6c5, modified to form only one star and use optimal 
SEDs. 
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at E = 2 keV. The second simulation employs the optimal 
four-bin SEDs found in Tables Q] and [2] 




(a) 




Figure 9. Comparison of hydrogen and helium ionization (a) and temper- 
ature profiles (b) around a power-law X-ray source after 50 Myr showing 
the solutions obtained using continuous (black) and optimal four-bin discrete 
(blue symbols) SEDs. 

We run two simulations, each identical to the other except 
for the choice of discrete SED. Our 'control' simulation uses 
monochromatic SEDs — the PopIII star is a monochromatic 
source of E = 29.6 eV photons, while the X-ray source emits 
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Figure 10. Comparison of the four-bin solutions of Wise & Abel] 1 1201 II) 
(black) and our own (blue crosses) in a radiation-hydrodynamic simulation 
using the Enzo code. The setup is the same as in RT06 Test Problem 2, ex- 
cept hydrodynamics is included. 



As shown in FigureQT| the magnitude of the errors between 
monochromatic and « v = 4 solutions is even more significant 
in the cosmological problem than in the RT06 test problem, 
since the ionizing luminosity of the blackbody source con- 
sidered is nearly two orders of magnitude larger (1 .2 x 10 50 
versus 5 x 10 48 s _1 ). For very luminous sources, even small 
errors in <t> and *P will become noticeable as characteristic 
timescales for photoionization and heating are short. 

During the BH phase of evolution, there are more ways 
for the monochromatic and multi-frequency solutions to differ 
aside from the SEDs being employed. The accretion luminos- 
ity depends on local gas properties, which will be different in 
each simulation due to errors accrued during the PopIII star's 
lifetime. Properties of the broader medium will of course vary 
for the same reason, leading to changes in how far soft X-rays 
are able to propagate before being absorbed. Throughout the 
100 Myr of evolution after the PopIII star's death, the Bondi- 
Hoyle accretion rate and thus luminosity of the accreting BH 
is on average an order of magnitude smaller in the nv = 4 
simulation than for the monochromatic case. Errors in ion- 
ization and temperature exceeding an order of magnitude per- 
sist throughout the BH phase as well. Rather than attempt 
to disentangle the BH phase induced errors from the preex- 
isting errors, we simply emphasize that SED-induced errors 
will compound in feedback situations like this, since the ini- 
tial conditions of each subsequent generation of objects will 
have been contaminated by errors associated with the previ- 
ous one. 

We cannot comment on the relative errors between 
monochromatic and multi-frequency treatments beyond the 
outermost column density contour, as our optimization ex- 
tended only to Nh i = 3.1 X 10 19 cm~ 2 . Future work focused 
on larger cosmological volumes, more luminous sources, and 
harder radiation fields will need to construct optimal SEDs 
valid beyond Nh i = 10 20 cm~ 2 , at least. 



6. DISCUSSION 

Algorithms developed for the purpose of studying point- 
source radiation (e.g., ray-tracing) are in principle capable 
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Figure 11. Ratio of slices of the ionized fraction (a) and temperature (b) ob- 
tained using our optimized n v = 4 blackbody SED (xh n 4 , T4) and the standard 
monochromatic SED (.vh iii , T\ ). Both slices are 2.25 Myr after the formation 
of a Population III star. Contours (from center outwards) correspond to hy- 
drogen column densities of Nn , = 2 and 4 x 10" cm -2 . 

of propagating continuous radiation fields, that is, tabulating 
Equations (TTOb and (fTTI) and computing ionization and heat- 
ing rates via Equations (fT2l i- (fl4l i. The reason many have not 
taken this approach could be due to the additional computa- 
tional overhead involved with using continuous SEDs — the 
quantities <t>,' and *F; must be tabulated over the complete col- 
umn density interval of interest. This includes column den- 
sities of all absorbing species, each of which must extend 
from the smallest expected column (i.e., the column density 
of a "fully ionized" cell — we adopted a minimum species 
fraction of jc m ; n = 10~ 5 ) up to the largest expected column 
(i.e., the column density of a fully neutral medium). The di- 
mensionality of 4>; and 1*, can be increased even further if 
for e xample energy-dependent secondary electron tr eatments 
(e.g.. lRicotti et al]|2002l: iFurlanetto & Stoeverj|2010l) or time- 
dependent SEDs are of interest. 

For the simulations of Section 13.21 we generated three- 
dimensional lookup tables for <!>, and covering the column 
density range 10 11 <N Hi < 10 21 , and 10 10 < Afaei^Heii < 
10 20 , sampling Nr i at 200 points, and A^e 1 and A^e 11 with 



100 points each, resulting in six three-dimensional tables, 
each consisting of 2 x 10 6 elements. We found that poorer 
sampling (e.g., tables of dimension 100 x 50 x 50) leads to 
artificial "notches" in ionization and temperature profiles due 
to errors in the bilinear interpolation. In our case, <I>h i = 
^He 1 = ^He n and 1 = ^He 1 = ^He 11 since all emission 
occurs above 10 2 eV, making the lower limit of integration 
for each quantity identical. In the general case, where emis- 
sion extends all the way to the hydrogen ionization threshold, 
all six quantities would be unique. Generating these tables 
can take hundreds of CPU hours or more for a single SED de- 
pending on the number of column density elements. In addi- 
tion, the radiative transfer solver requires additional modules 
to read in the lookup table, and perform interpolation four 
times per absorbing species per grid element (see Eqs (fl2l - 
(fT4l>). For sources with discrete SEDs, one can simply com- 
pute the photo-ionization rate for each neutral species, from 
which point the secondary ionization and heating rate coef- 
ficients are obtained in a simple algebraic fashion (see Eqs 

G3J-Q3). 

For high-resolutio n simulations focu s ed on a single source 
of ra diation (e.g., iKuhlen & Madaul l2005t lAlvarez et ail 
2009), the additional effort required to accommodate con- 
tinuous radiation fields seems well worth it to ensure that 
the ionization and thermal state of the gas is captured accu- 
rately. However, in large-scale simulations of cosmic reion- 
ization, which may spawn hundreds of thousands or perhaps 
millions of radiating 'star particles' (depending on the simu- 
lation volume, resolution, etc.), ray-tracing methods are cer- 
tainly not the most computationally advantageous algorithm. 
This is because the computational cost of a ray-tracing cal- 
culation scales with the number of radiation sources and the 
number of frequency bins in each source SED (though the 
former c ost can be mitiga t ed by merging nearby radiation 
sources; iTrac&Cenl 1200 7; Oka moto et all I2012D . If pho- 
tons with long mean free paths are of interest, the simula- 
tion will be even more expensive since rays must be fol- 
lowed to larger distances, i.e., more ray segments and it- 
erations of the numerical solver are required. An appeal- 
ing option is to instead use moment-based methods such as 
the Va riable Eddington Ten sor approach (e.g jGnedin & Abel 
200UlPetkova & Srjringei l2009l) . flux-limited diffusion (e.g. 
Reynolds et al. 2009)," o r other variations (Gonzal ez et al.l 
2007t lAubert & Teyssieri 12008b iFinlator et al l 120091) . as the 
computational cost of such algorithms is independent of the 
number of radiation sources and the mean free paths of pho- 
tons, scaling only with the number of frequency bins in each 
source spectrum. 

As discussed in SectionQ] multi-group schemes common in 
the literature are an improvement over fiducial discrete SEDs, 
though it is not generally clear how many bandpasses are re- 
quired for a given problem, or where they should lie in fre- 
quency space. Moreover, multi-group radiation suffers from 
the same problem as discrete polychromatic emission: pho- 
tons at each frequency are absorbed near a characteristic col- 
umn density, A^har- Computing new spectrum-weighted ab- 
sorption cross-sections, a„, for each frequency group merely 
shifts the location of Ach ar . 

In principle, our minimization technique could be used to 
optimally select which bandpasses should be used for a multi- 
group algorithm, though in practice it would be much more 
computationally expensive. Rather than varying the location 
(v„) or normalization (/„) of frequency bin n on each Monte 
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Carlo step, one would instead vary the position of bandpass 
edges, which would change the mean photon energy in each 
bandpass (hv„) and spect rum-weighted cross section, d„ (e.g., 
lAubert & Teyssierl 120081) . Because hv„ and <7„ are integral 
quantities, they would need to be computed numerically on 
each Monte-Carlo step, and thus hundreds of thousands of 
times for a single optimization. 

7. CONCLUSIONS 

We have shown that the manner in which a discrete SED 
is constructed can induce substantial errors in simulation re- 
sults, both in the ionization and temperature profiles around 
stars and quasars. But, these errors can be avoided to a large 
degree using only four discrete emission frequencies if source 
SEDs are designed via the methods of Section [4] Discrete 
SEDs constructed in a simple way (e.g., bins linearly spaced 
in frequency) will perform more poorly than optimally se- 
lected SEDs with the same number of bins, since it is the 
column density interval of interest that dictates the range of 
photon energies required, and the power to which each is as- 
signed. 

In general, discrete SED treatments fail to ionize and/or 
heat gas at large column densities, i.e., large physical scales 
or environments with dense clumps of gas. This has strong 
implications for simulations dedicated to understanding the 
magnitude and mode of radiative feedback on gas surround- 
ing radiation sources. Current questions of this sort include 
whether or not radiation stimulates or suppresses further star 
formation in nearby proto-stellar clouds, and if radiative feed- 
back can stifle the growth of SMBHs at high redshift. 

As expected, extending our one-dimensional work to three- 
dimensions produces ionized regions around a first star and 
remnant BH that deviate significantly in ionized fraction, tem- 
perature, size, and morphology. Such findings have implica- 
tions in radiative feedback, but also in studies of both hydro- 
gen and helium reionization. Certainly miscalculations of the 
ionization state of gas surrounding galaxies in the early uni- 
verse will lead to errors in the volume averaged neutral frac- 
tion, volume filling factor of ionized gas, and the optical depth 
of the CMB to electron scattering (x e ). As we demonstrated 
in Section[3] such errors also introduce uncertainties in the in- 
terpretation of future 21 cm measurements, since the primary 
observable quantity (87i) depends directly on the hydrogen 
neutral fraction, electron density, and gas kinetic temperature. 

Our optimizations in this work are by no means com- 
prehensive, having selected two commonly used radiation 
sources (UV blackbody and X-ray power law) as test cases 
to demonstrate the method. However, optimization for more 
complex spectra is straightforward, and any new optimiza- 
tions run will be made publicly available by the authors. The 
minimization code and one-dimensional radiative transfer 
codes are both available upon request. We leave more 
detailed investigations of reionization and radiative feedback, 
including multiple radiation sources and multi-frequency 
radiation transport, to future work. 

The authors thank Steven Furlanetto and Daniel Reynolds 
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as the anonymous referee for a thorough review and 
many helpful suggestions. The LUNAR consortium 
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APPENDIX 

A. OPTIMIZATION VIA SIMULATED ANNEALING 

To solve Equation (I24K we employ the Monte Carlo method of Simulated Annealing (Kirkpa trick et alJll983t ICerny||1985l) . 
For a given source and n v , we run K Monte-Carlo trials, each consisting of L steps, aimed at determining the optimal values of 
/„ and v„ for n v frequency bins. We do not require the bolometric luminosity of sources to be conserved (i.e., Yl^ = \hi 7^ 1 is 
allowed), since some photons may traverse the entire one-dimensional "volume" without ionizing a single atom, or some fraction 
of the luminosity may be emitted below the hydrogen ionization threshold. Inclusion of such photons would be computational 
effort wasted in a fully three-dimensional ray-tracing calculation, for example, since their mean free paths are very long, and 
once absorbed they may contribute negligibly to ionization and heating. 

Each random walk begins with randomly generated values of v„ distributed between the hydrogen ionization threshold and the 
maximum emission frequency in the spectrum, and randomly generated values of /„ that sum to unity. Subsequent steps vary the 
energy or normalization of (randomly chosen) frequency bin n. In order to steer each random walk towards the global minimum, 
we first evaluate the quantity 

P = exp[-(/jy -/*,/_! )/7sa] (Al) 

where k = 0, 1,2, . . . ,K represents the current step in the current random walk, /, where I = 0, 1,2,... ,L, and / is the "cost 
function," a measure of how good our current solution is. We adopt a cost function which is the sum of errors in <!>,■ and *?, over 
the column density range of interest. For each species (/), and each integral quantity (<f>, V P), we add the maximum deviation from 
continuous and discrete solutions in the opticall y thi n limit (first term in Equation (IA2t ). the maximum deviation over the entire 
column dens ity r ange (second term in Equation (IA2t ). and the average deviation over the entire column density range (final term 
in Equation (IA2I) ). all in dex, i.e., 



i A=4>. l P 



max 



log 



A, 



Aj(v^/,4,/) 



x=0- 



log 



A^Vjt,/,/*,/) 



x>0. 



A^Vjt,/,/*,/) 



T>0 



At each step in a given random walk, we also generate a random number, q £ [0,1], that will determine whether we keep our 
current guess, {Vk.hh.i), or revert to our previous guess, (v*,/_i ,4./-i )■ The condition for keeping our current guess is P > q. 

The key aspect of this analysis is how we vary the control parameter 7sa, which is called the temperature in analogy with 
Boltzmann's equation (we add the s ubscript SA to distinguish the gas kinetic temperature from this unphysical Simulated An- 
nealing temperature). Equation (IA1 b tells us that regardless of the value of Tsa, if fk.i < fk,i-i (i.e., our most recent guess is better 
than the last), then P > 1, and we have a 100% chance of keeping our current guess. In other words, our method of controlling 
the Tsa only effects how we deal wi th bad guesses — decreasing the temperature means we become less tolerant of bad guesses. 
There are many ways of doing this (Pre ss et ailll992l) . but for simplicity we adopt the following technique. Every s/n v steps per 
frequency bin, we take 

T^XT, (A3) 

where X is an experimentally determined quantity of order unity. For all results presented here, we have adopted X = 0.98, and 
s/riy = 10. We change the number of steps per random walk depending on the dimensionality, 2n v . We have found through 
experimentation that a good rule of thumb is L = 5000 steps per trial, K, per frequency bin n v for our choice of X and s /n v . These 
control parameters are fairly conservative — further experimentation with them may yield converged solutions for fewer trials, 
K, and steps, L. 

B. CODE VERIFICATION 

Our one-dimensional radiative transfer code solves Equations <[TJ — <SJ using the implicit Euler method for integration and a 
Newton-Raphson technique for root finding. Each simulation is initialized on a grid of N c cells between Lq and Lb m , such that 
the finest resolution element is Ax = (L(,ox — Lo)/N c , or simply Ax = 1 /N c in code units. Gas inside of the start radius, Lq, 
contributes no optical depth, and Equations (Q]i-© are not solved. For the purposes of this section, we chose to use N c linearly 
spaced cells between Lq and Lbox, though our code allows arbitrarily structured grids. 

In order to track the propagation of ioniza tion fronts accurately, we limit the time-step based on a maximum neutral fraction 
change as introduced in Shapi ro et al.1 d2004), 

Ar,=e ion ' -7 , (Bl) 
\ani/at\ 

where we include all absorbing species , ; =H I, He I , He II , and set At = min(Af,). We additionally require that the time step 
increase by a factor of two at most, as in Wis e & Abell (1201 lh . For all simulations presented in this work, we have set £i on = 0.05. 

The primary solver implemented in our code assumes the speed -of-light is infinite. Such an algorithm is appealing for two 
main reasons, aside from the fact that it is a very good approximation for the problems presented in this work. First, treating the 
speed-of-light explicitly introduces additional computational overhead as "photon packages" must be launched from the radiation 
source at each time step and tracked unti l they exit the domain. In the earliest stages of I-front propagation, the time step can be 
very small (as required by Equation (IBU ). meaning the total number of photon packages, N p , will be much larger than the total 
number of grid cells, N c . Whereas c = °° treatments only require Equations ([T}-(|4) to be solved once per cell, finite speed-of-light 
treatments require this system of equations to be solved for each photon package. At later times, when N p < N c , solving the 
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ion and heat equations is cheaper for finite speed-of-light treatments, though this offers no real advantage since the majority 
of the computational expense is at early times when I-front propagation is fastest. We have also included a finite c solver to 
accommodate a broader class of problems that may be of interest in future work. 

The second advantage of assuming c = °° is that it allows the code to be efficiently parallelized. If c — °°, cells in the domain 
can be solved in arbitrary order by a single processor, or simultaneously by a network of processors, since the radiation incident 
on any cell is predetermined at the outset of each individual time step. Previous authors have ensured causality by solving cell k 
before cell lc+ 1 at time t (where increasing k corresponds to increasing r), but this is not in fact necessary — causality is ensured 
by the monotonicity of column density with distance. In other words, when c = °°, Ni does not change within any given time 
step, and so the column density (and thus radiative flux) to cell k is less than the column density (and flux) to cell k+1, meaning 
the solution of Equations (HJ-© in cell k + 1 is completely independent of the properties of cell k at time t + At . 

To d emonstrate the functionality of the code, we repeat tests 1 and 2 from the Radiative Transfer Comparison Project (llliev et alJ 
( 2006, hereafter referred to as RT06)) on a grid of 200 linearly spaced cells. Test 1 is the expansion of an H II region in a hydrogen- 
only, isothermal medium surrounding a monochromatic source of 13.6 eV photons. We adopt the same parameters used in RT06: 
constant temperature T = 10 4 K, uniform hydrogen number density «h = 10~ 3 cm~ 3 , ionized fraction *h n = 1.2 X 10 3, in a 
box Lbox = 6.6 kpc in size, and with photon luminosity Q = 5 x 10 48 s . The classical analytic solution for the radius of an 
ionization front is 

r !F (t) = r s (l-e- t ^) 1 / 3 , (B2) 

where r s is the Stromgren radius, 

( 36 

and the recombination time, f rec , is defined as 

f rec = • (B4) 

O^H ii«H i 

This solution is approximate eve n in isothermal media, given that it assumes a constant neutral hydrogen density, «h i- More 
accurate analytic solutions exist (lOsterbrock & Ferlandll2006T) . and predict a departu re from the classical solution at f/f rec — 1, 
which grows to a ~ 5% difference by f/f rec — 4. Our numerical solution (see Figure [12(a)} captures this behavior very well. In 
Figure [T2(b)| we show radial profiles of the ionized and neutral fractions at three stages of the I-front expansion, which are again 
in very good agreement with the calculations presented in RT06. 





Figure 12. Test 1 : (a) Comparison of the numerical (dashed) and analytic (solid) solutions for the position of an expanding ionization front as a function of time 
in a hydrogen-only, isothermal medium (RT06 problem 1; top), and the ratio of the calculated and analytic solutions as a function of time and grid resolution 
(bottom). The numerical solution displayed in the top panel is from the highest resolution simulation (800 grid cells, i.e., Ax = L[, ox /800). ((b)) Radial profiles 
of the neutral (solid) and ionized (dashed) fractions at t = 10, 100, and 500 Myr. 

Test 2 is the same as Test 1, except now the temperature is allowed to evolve according to Equation |@}, and the monochromatic 
radiation source is replaced by a 10 5 K blackbody spectrum. Radial profiles of the neutral and ionized fractions and temperature 
can be seen in FigureQj] Again, our numerical solutions are in very good agreement with previous work. 
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